Nucleon form factors and moments of generalized parton distributions using 

Nf = 2 + 1 + 1 twisted mass fermions 
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We present results on the axial and the electromagnetic form factors of the nucleon, as well as, 
on the first moments of the nucleon generalized parton distributions using maximally twisted mass 
fermions. We analyze two Nf— 2+1+1 ensembles having pion masses of 10 MeV and 354 MeV 
at two values of the lattice spacing. The lattice scale is determined using the nucleon mass com- 
puted on a total of 18 Nf =2+1+1 ensembles generated at three values of the lattice spacing, a. 
The renormalization constants are evaluated non-perturbatively with a perturbative subtraction of 
0(a 2 )-terms. The moments of the generalized parton distributions are given in the MS scheme at a 
scale of = 2 GeV. We compare with recent results obtained using different discretization schemes. 
The implications on the spin content of the nucleon are also discussed. 

PACS numbers: 11.15.Ha, 12.38.Gc, 12.38.Aw, 12.38.-t, 14.70.Dj 
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I. INTRODUCTION 

Recent progress in the numerical simulation of Lattice 
Quantum Chromodynamics (LQCD) has been remark- 
able. The improvements in the algorithms used and the 
increase in computational power have enabled simula- 
tions to be carried out at near physical parameters of the 
theory. This opens up exciting possibilities for ab ini- 
tio calculation of experimentally measured quantities, as 
well as, for predicting quantities that are not easily ac- 
cessible to experiment. Understanding nucleon structure 
from first principles is considered a milestone of hadronic 
physics and a rich experimental program has been de- 
voted to its study, starting with the measurements of 
the electromagnetic form factors initiated more than 50 
years ago. Reproducing these key observables within the 
LQCD formulation is a prerequisite to obtaining reliable 
predictions on observables that explore Physics beyond 
the standard model. 

A number of major collaborations have been studying 
nucleon structure within LQCD for many years. How- 
ever, it is only recently that these quantities can be ob- 
tained with near physical parameters both in terms of the 
value of the pion mass, as well as, with respect to the con- 
tinuum limit (lMTlj . The nucleon electromagnetic form 
factors are a well suited experimental probe for study- 
ing nucleon structure and thus provide a valuable bench- 
mark for LQCD. The nucleon form factors connected to 
the axial- vector current are more difficult to measure and 
therefore less accurately known than its electromagnetic 



form factors. A notable exception is the nucleon axial 
charge, gA, which is accurately measured in /3-decays. 
The fact that gA can be extracted at zero momentum 
transfer and that it is technically straight forward to be 
computed in LQCD, due to its isovector nature, makes it 
an ideal benchmark quantity for LQCD. The Generalized 
Parton Distributions (GPDs) encode information related 
to nucleon structure that complements the information 
extracted from form factors [12Hl4l |. They enter in sev- 
eral physical processes such as Deeply Virtual Compton 
Scattering and Deeply Virtual Meson Production. Their 
forward limit coincides with the usual parton distribu- 
tions and, using Ji's sum rule [HI, allows one to deter- 
mine the contribution of a specific parton to the nucleon 
spin. In the context of the "proton spin puzzle" , which 
refers to the unexpectedly small fraction of the total spin 
of the nucleon carried by quarks, this has triggered in- 
tense experimental activity |16l - [20l |. 



II. LATTICE EVALUATION 

In this work we consider the nucleon matrix elements 
of the vector and axial- vector operators 
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FIG. 1: Connected nucleon three-point function. 



where r a are the Pauli matrices acting in flavor space, 
ip denotes the quark doublet with components the up- 
and down-quarks. In this work we consider the isovector 
combination by taking a — 3, except when we discuss 
the spin fraction carried by each quark. Furthermore, we 
limit ourselves to n = 1 and n — 2. The case n = 1 
reduces to the nucleon form factors of the vector and 
axial-vector currents, while n — 2 correspond to matrix 
elements of operators with a single derivative. The curly 
brackets represent a symmetrization over indices and sub- 
traction of traces, only applicable to the operators with 
derivatives. There are well developed methods to com- 
pute the so called connected diagram, depicted in Fig. [TJ 
contributing to the matrix elements of these operators 
in LQCD. Each operator can be decomposed in terms of 
generalized form factors (GFFs) as follows: The matrix 
element of the local vector current, Oy 3 , is expressed as 
a function of the Dirac and Pauli form factors 



(N(p',s')\0» s \N{p,s)) = 
un(p',s') 
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where un{p, s) denote the nucleon spinors of a given mo- 
mentum p and spin s. fi(0) measures the nucleon charge 
while F 2 (0) measures the anomalous magnetic moment. 
They are connected to the electric, Ge, and magnetic, 
Gm, Sachs form factors by the relations 
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The local axial current matrix element of the nucleon 
(N(p', s')|0^3 1 N(p, s)) can be expressed in terms of the 
form factors G a and G p as 



(N(j/,a')\CK a \N(p,a)) 



u N (p',s' 



~u N (p,s) .(4) 



The matrix elements of the one derivative op- 
erators are parameterized in terms of the GFFs 
A20 {q 2 ), B 20 (q 2 ), C 20 (q 2 ), and A 20 (q 2 ) and B 20 {q 2 ) for 
the vector and axial-vector operators respectively, ac- 
cording to 



(N(p',s')\0^\N(p,s)) = u N (p', s') A 20 (q 2 ) 7 { ^ } + B 20 (q 2 ) + c 20 ( q 2 ) -q^] \u N {p, s) , (5) 

771 J z 



2m 

(N(p>, s')\0% \N(p, s)) = u N (p', s') \A 20 (q 2 ) 7 ^F-> 7 5 + B 20 (q 2 ) ^^7 5 1 \u N {p, s) . 

L 2m J 2 



(6) 



Note that the GFFs depend only on the momentum 
transfer squared, q 2 — (p' — p) 2 , p' is the final and p 
the initial momentum. The isospin limit corresponds 
to taking r 3 /2 in Eq. © and gives the form factor 
of the proton minus the form factors of the neutron. 
In the forward limit we thus have Ge(0) = 1 and 
Gm(0) = [i v — [i n — \ = 3.7 [2l[, which is the isovec- 
tor anomalous magnetic moment. Similarly, we obtain 
the nucleon axial charge, Ga(0) = <?a, the isovector mo- 
mentum fraction, A 2 q (0) = (x) u -d and the moment of 
the polarized quark distribution, ^20(0) = (x)au— Ad- I n 
order, to find the spin and angular momentum carried 
by each quark individually in the nucleon we need the 
isoscalar axial charge and the isoscalar one-derivative ma- 



trix elements of the vector operator. Unlike the isovector 
combinations, where disconnected fermion loops vanish 
in the continuum limit, the isoscalar cases receive con- 
tributions from disconnected fermion loops. These con- 
tributions are believed to be small and are usually ne- 
glected. Their evaluation is difficult due to the compu- 
tational cost but techniques are being developed to com- 
pute them [23| . If one assumes that the disconnected 
contributions are small, it is straightforward to evaluate 
the isoscalar matrix elements taking into account only 
the connected part depicted in Fig. [TJ The quark con- 
tribution to the nucleon spin is obtained using Ji's sum 
rule: J q = |[A| (0) + -Bf (0)]. Moreover, using the axial 
charge for each quark, g q Al we obtain the intrinsic spin 
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of each quark, AE 9 = g q A , and via the decomposition 
J q = + L q we can extract the quark orbital angu- 

lar momentum L q . 



In the present work we employ the twisted mass 
fermion (TMF) action [23| and the Iwasaki improved 
gauge action [2J]. Twisted mass fermions provide an at- 
tractive formulation of lattice QCD that allows for auto- 
matic 0(a) improvement, infrared regularization of small 
eigenvalues and fast dynamical simulations (25[. In the 
computation of GFFs the automatic 0{a) improvement 
is particularly relevant since it is achieved by tuning only 
one parameter in the action, requiring no further im- 
provements on the operator level. 

We use the twisted mass Wilson action for the light 
doublet of quarks 

Si = ^Xl( x ) [ D w+m {0t i)+ij 5 T 3 fit] xi(x) , (7) 

X 

where Dw is the Wilson Dirac operator, W(o,i) is the un- 
twisted bare quark mass, fii is the bare light twisted mass. 
The quark fields xi are m the so-called "twisted basis" 
obtained from the "physical basis" at maximal twist by 
the transformation 

ip=-^=[l + iT 3 j 5 ] X i and iP=xi-^=[l + iT 3 l5 }. (8) 
In addition to the light sector, we introduce a twisted 



heavy mass-split doublet \h = (Xc,Xs) for the strange 
and charm quarks, described by the action 



Sh = ^Xhix) [ D w+m{o,h)+il<iT 1 fL a + T 3 fi 5 ] Xh(x) , 

X 

(9) 

where m^ o h ^ is the untwisted bare quark mass for the 
heavy doublet, /i CT is the bare twisted mass along the r 1 
direction and fis is the mass splitting in the r 3 direction. 
The quark mass rn(o,/i) is set equal to Ti( 0j j) in the simula- 
tions thus ensuring C(a)-improvement also in the heavy 
quark sector. The chiral rotation for the heavy quarks 
from the twisted to the physical basis is 
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[1 + ir j 5 ]xh and ip=x h -^ [l+ir 1 ^]. (10) 



The reader can find more details on the twisted mass 
fermion action in Ref. 26] . 



A. Correlation functions 

The GFFs are extracted from dimensionless ratios 
of correlation functions, involving two-point and three- 
point functions that are defined by 
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For the insertion, O^ 1 '"^", we employ the vector 
(-07^), the axial-vector (05 7 s 7^0), the one-derivative 
vector (0 ^f^ 1 ip) and the one-derivative axial- vector 
(-0 75 7^ 1 Z? M2 '' V) operators. We consider kinematics for 
which the final momentum p 1 = and in our approach we 
fix the time separation between sink and source tf — U. 
The projection matrices T° and T k are given by 



r° = ~(i + -*>), 



1 «757fe • 



(13) 



The proton interpolating field written in terms of the 
quark fields in the twisted basis (u and d) at maximal 
twist is given by 



J{x) 



1 



a be 



u aT > 



(z)C 75 (?(x) u c {x), (14) 



where C is the charge conjugation matrix. We use Gaus- 
sian smeared quark fields [23, to increase the overlap 
with the proton state and decrease overlap with excited 
states. The smeared interpolating fields are given by 



(15) 
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H(x,y;U(t)) 
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We also apply APE-smearing to the gauge fields enter- 
ing the hopping matrix H. The parameters for the Gaus- 
sian smearing an and Na are optimized using the nucleon 
ground state [29( . Different combination of Gaussian pa- 
rameters, Nq and og, have been tested and it was found 
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that combinations of Nq and olg that give a root mean 
square radius of about 0.5 fm are optimal for suppress- 
ing excited states. The results of this work have been 
produced with 

P = 1.95 : N G = 50 , a G = 4, A APE = 20, &ape = 0.5, 
P = 2.10 : N G = 110, a G = 4, TVape = 50, &ape = 0.5 . 

As already point out, in correlators of isovector operators 
the disconnected diagrams are zero up to lattice artifacts, 
and can be safely neglected as we approach the contin- 
uum limit. Thus, these correlators can be calculated by 
evaluating the connected diagram of Fig. [T] for which 
we employ sequential inversions through the sink [30| . 
The creation operator is taken at a fixed position Xi=0 
(source). The annihilation operator at a later time tf 
(sink) carries momentum p 1 =0. The current couples to a 
quark at an intermediate time t and carries momentum 
q. Translation invariance enforces q = —p for our kine- 
matics. At a fixed sink-source time separation we obtain 
results for all possible momentum transfers and insertion 
times as well as for any operator Q^ 1 '"^ ) with one 
set of sequential inversions per choice of the sink. We 
perform separate inversions for the two projection ma- 
trices r° and J2k^ k gi ven m Eq. (fT5|) . An alternative 
approach that computes the spatial all-to-all propagator 
stochastically has shown ot be suitable for the evalua- 
tion of nucleon three-point functions [3l|. Within this 
approach one can include any projection without need- 
ing additional inversions. 

Using the two- and three-point functions of Eqs. (fTTj) - 
(|12p and considering operators with up to one derivative 
we form the ratio 



R^{T x ,q,t) = 



G^(T x ,q,t) 
G0,t f -ti) 



x G(p,t f -t)G(0,t-U)G(0,t f -U) 
V G(6,t f -t)G(f,t-U)G(f,t f -U) A 

which is optimized because it does not contain poten- 
tially noisy two-point functions at large separations and 
because correlations between its different factors reduce 
the statistical noise. For sufficiently large separations 
tf — t and t — ti this ratio becomes time-independent 
(plateau region): 



lim 

— t— >oc t- 



lim 

t;-K 



R^(T x ,q,t) =U^(r x ,q). (17) 



From the plateau values of the renormalized asymptotic 
ratio H(r x ,q)fi = ZH(T x ,q) the nucleon matrix ele- 
ments of all our operators can be extracted. The equa- 
tions relating tt(T x ,q) to the GFFs can be found in 
Refs. |1"3]. All values of q resulting in the same q 2 , the 
two choices of projector matrices r° and J^k ^ k gi ven 
Eq. (|13l) and the relevant orientations p., v of the oper- 
ators lead to an over-constrained system of equations, 
which is solved in the least-squares sense via a singular 



value decomposition of the coefficient matrix. All quan- 
tities will be given in Euclidean space with Q 2 = —q 2 the 
Euclidean momentum transfer squared. Both projectors 
r° and J2k are required to obtain all GFFs, except 
for the case of the local axial-vector operator, for which 
the projection with T° leads to zero. For the one deriva- 
tive vector operator, both cases /i = v and ji ^ v are 
necessary to extract all three GFFs, which on a lattice 
renormalizc differently from each other [32| . On the other 
hand, the one-derivative axial-vector form factors can be 
extracted using only correlators with /i ^ v, but we use 
all combinations of p, v in order to increase statistics. In 
Fig. [5] we show representative plateaus for the ratios of 
the local axial- vector and the one derivative vector oper- 
ators at P = 1.95, using different momenta, projectors, 
and indices fi, v. 




-0.15 - 
0.05 — 

0.04 - 

0.03 - 



■ ■ ■ ■ ■ 



-0.05 - 
-0.1 - 
-0.15 - 



g 0.04 

£ 0.03 

4 0.02 

%. 0.01 L 



FIG. 2: Ratios for the matrix elements of the local axial- 
vector operator (upper) and one derivative vector operator 
(lower) for a few exemplary choices of the momentum. The 
solid lines with the bands indicate the fitted plateau values 
with their jackknife errors. From top to bottom the momen- 
tum takes values p = (0,0,0), (1,0,0), (0,-1,0) and (1,0,1). 



Since we use sequential inversions through the sink we 
need to fix the sink-source separation. Optimally, one 
wants to keep the statistical errors on the ratio of Eq. (fTBj) 
as small as possible by using the smallest value for the 
sink-source time separation that still ensures that the ex- 
cited state contributions are sufficiently suppressed. Re- 
cent studies have shown that the optimal sink-source sep- 
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aration is operator dependent [33j, |34[. For g& excited 
state contamination was found to be small. We have 
also tested different values of the sink-source time sepa- 
ration [3J] for the magnetic form factor and found consis- 
tent results when the sink-source separation was about 
1 fm within our statistical accuracy. For the momentum 
fraction one would need to re-examine the optimal sink- 
source separation, which would require a dedicate high 
accuracy study. Since in this work we are computing 
several observables, we will we use tf — ti ~ 1 fm that 
correspond to the following values 

,9 = 1.95: (t f -ti)/a=12, ft = 2.10 : (t f - t t )/a=18. 

This choice allows to compare with other lattice QCD 
results where similar values were used. 



B. Simulation details 

In Table U we tabulate the input parameters of the cal- 
culation, namely ft, L/a and the light quark mass aa, as 
well as, the value of the pion mass in lattice units |35ll36j. 
The strange and charm quark masses were fixed to ap- 
proximately reproduce the physical kaon and D-meson 
masses, respectively [37j- The lattice spacing a given in 
this Table is determined from the nucleon mass as ex- 
plained in the following subsection and it will be used 
for the baryon observables discussed in this paper. We 
note that the study of the systematic error in the scale 
setting using the pion decay constant as compared to the 
value extracted using the nucleon mass is currently be- 
ing pursued. Since the GFFs are dimensionless they are 
not affected by the scale setting. However, a is needed 
to convert Q 2 to physical units and therefore it does af- 
fect quantities like the anomalous magnetic moment and 
Dirac and Pauli radii since these are dimensionful pa- 
rameters and depend on fitting the Q 2 -dependence of the 
form factors. 



ft = 1.95, a = 0.0863(11) fm, r Q /a = 5.66(3) 



32 3 x 64, L = 2.5 fm 


an 


0.0055 




No. of confs 


950 




am.Tr 


0.15518(21)(33) 




LlTlr: 


4.97 



ft = 2.10, a = 0.0657(7) fm, ro/a = 7.61(6) 



48 3 x 96, L = 2.9 fm 




0.0015 




No. of confs 


900 




a v/Itx 


0.06975(20) 






3.35 



TABLE I: Input parameters (ft, L, a/x) of our lattice calcu- 
lation with the corresponding lattice spacing a, determined 
from the nucleon mass, and pion mass am^ in lattice units. 



C. Determination of lattice spacing 

For the observables discussed in this work the nucleon 
mass at the physical point is the most appropriate quan- 
tity to set the scale. The values for the nucleon mass 
were computed using 7V/=2+l+l ensembles for /3=1.90, 
/3=1.95 and /3=2.10, a range of pion masses and volumes. 
To extract the mass we consider the two-point correlators 
defined in Eq. (fTTj) and construct the effective mass 
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where Aj = Ej 



rriN is the energy difference of the 
excited state j with respect to the ground state mass, 
toat. Our fitting procedure to extract is as follows: 
The mass is obtained from a constant fit to mf(t) for 
t > t\ for which the contamination of excited states is 
believed to be small. We denote the value extracted as 



i (A) 

N 



(ti). A second fit to m e ^(t) is performed including 
the first excited state for t > t\, where t[ is taken to be 
2a or 3a. We denote the value for the ground state mass 
extracted from the fit to two exponentials m N . We vary 
t\ such that the ratio 
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where 

N 

(A) u s . (B) 

am x N (ti) + am x N 



(19) 



drops below 50% of the statistical error on m^(ti). The 
resulting values for the nucleon mass are collected in Ta- 
ble HI 

In Fig. [3] we show results at three values of the lattice 
spacing corresponding to ft=1.90, /3=1.95 and /3=2.10. 
As can be seen, cut-off effects are negligible and we can 
therefore use continuum chiral perturbation theory to ex- 
trapolate to the physical point using all the lattice re- 
sults. 

To chirally extrapolate we use the well-established 
0(p 3 ) result of chiral perturbation theory (xPT) given 
by 



m N 



4ci m„ 



3gj , 



(20) 



We perform a fit to the results at the three ft values given 
in Table |ll]and extract 



a^i.go = 0.0950(11)(20), 
00=1.98 = 0.0863(9)(20), 
ap=2.w = 0.0657(7) (15). 



(21) 



We would like to point out that our lattice results show 
a curvature supporting the m^-term. In order to esti- 
mate the systematic error due to the chiral extrapola- 
tion we also perform a fit using heavy baryon (HB) %PT 
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TABLE II: Values of the nucleon mass and the associated 
statistical error. 



to 0(p 4 ) with explicit A degrees of freedom in the so 
called small scale expansion(SSE) We take the dif- 
ference between the 0(p 3 ) and 0(p 4 ) mean values as an 
estimate of the uncertainty due to the chiral extrapo- 
lation. This error is given in the second parenthesis in 
Eqs. (f2"Tj) and it is about twice the statistical error. The 
values of the lattice spacing given in Eqs. (l2~Tj) will be 
used for converting to physical units the quantities we 
study here. We would like to point out that redoing the 
0(p 3 ) fit eliminating data for which Lm^ < 3.5 yields 
£10=1.90 = 0.0942(14) fm, a^ =1 . 95 = 0.08 58(11) fm and 
0-13=2.1 — 0.0653(8), which are consistent with the values 
given in Eq. (|21[) . In performing these fits we only take 
into account statistical errors. Systematic errors due to 
the choice of the plateau are not included. We also note 
that the lattice spacings were also determined from the 
pion decay constant using NLO SU(2) chiral perturba- 
tion theory to extrapolate the lattice data. The values 
obtained at (3 = 1.90, 1.95 and 2.10 in this preliminary 
analysis that included only a subset of the ensembles used 
here are smaller [3(| , as compared to the values extracted 
using the nucleon mass. For the two (3- values studied in 
this work they were found to be a/, = 0.0779(4) fm at 
(3 = 1.95 and a/„ = 0.0607(3) fm at (3 = 2.10, where 
with a/ T we denote the lattice spacing determined us- 
ing the pion decay constant. This means that the values 
of the pion mass in physical units quoted in this paper 
are equivalently smaller than those obtained using a to 
convert to physical units. A comprehensive analysis of 
the scale setting and the associated systematic uncertain- 
ties is currently being carried out by European Twisted 
Mass Collaboration (ETMC) and will appear elsewhere. 
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FIG. 3: Nucleon mass at three lattice spacings. The solid lines 
are fits to 0(p 3 ) and 0(p 4 ) HBxPT with explicit A degrees 
of freedom in the so called small scale expansion(SSE). The 
dotted lines denote the error band. The physical point is 
shown with the asterisk. 



D. Renormalization 

We determine the renormalization constants needed for 
the operators discussed in this work in the RI'-MOM 
scheme I38II by employing a momentum source at the 
vertex [39jj . The advantage of this method is the high 
statistical accuracy and the evaluation of the vertex for 
any operator including extended operators at no signifi- 
cant additional computational cost. For the details of the 
non-perturbative renormalization see Ref. [38] . In the RI 
scheme the renormalization constants are defined in the 
chiral limit. Since the mass of the strange and charm 
quarks are fixed to their physical values in these simu- 
lations, extrapolation to the chiral limit is not possible. 
Therefore, in order to compute the renormalization con- 
stants needed to obtain physical observables, ETMC has 
generated Nf=A ensembles for the same (3 values so that 
the chiral limit can be taken [4(|. Although we will use 
the Nf=4 ensembles for the final determination of the 
renormalization constants, it is also interesting to com- 
pute the renormalization constants using the Nf —2+1+1 
ensembles and study their quark mass dependence. This 
test was performed on both the (3 = 1.95 and the f3 = 2.10 
ensembles. In the upper panel of Fig. 2] we show results 
at P = 2.10 for both Nf=4 and iV/=2+l+l ensembles 
for the one derivative Z-factors in the RI'-MOM scheme. 
As can be seen, we obtain compatible values for all four 
cases. We also observe the same agreement for Zy and 
Za also at (3 — 1.95. This can be understood by examin- 
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FIG. 4: Upper panel: One derivative renormalization func- 
tions for /3 — 2.10, afi — 0.0015 using Nf=4 gauge configura- 
tions, where Zdvi (Zdai) = Z D y (Z D \) and Zdv2 (Zda2) = 
Z^y" (Z DA V ). Black circles are the unsubtracted data and 
the magenta diamonds the data after subtracting the pertur- 
bative C(a 2 )-terms. For comparison, we show the subtracted 
data using 7V/=2+l+l gauge configurations at the same value 
of the quark mass and /3 (blue crosses). Lower panel: One 
derivative renormalization functions for /3 = 1.95 using iV/=4 
gauge configurations as a function of the twisted quark mass. 



ing the quark mass dependence of these renormalization 
constants. In the lower panel of Fig. @]we show, for the 
Nf=4 case, the dependence of Zdv, Zda on four light 
quark masses. The values we find are consistent with 
each other. This explains the fact that the results in 
the Nf=4 and Nf =2+1+1 cases arc compatible. Fur- 
thermore, it makes any extrapolation of iV/=4 results to 
the chiral limit straight forward. We p erform a per- 
turbative subtraction of 0(a 2 )-terms [11 El, HIJ. This 
subtracts the leading cut-off effects yielding, in general, a 
weak dependence of the renormalization factors on (ap) 2 
for which the (ap) 2 -+ limit can be reliably taken, as 
can be seen in Figs. 0] and [5] for the two Nf = 2 + 1 + 1 
ensembles. We also take the chiral limit, although the 
quark mass dependence is negligible for the aforemen- 
tioned operators. 
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data are shown in black circles and the data after the 0(a 2 )- 
terms have been subtracted are shown in magenta diamonds. 
The solid diamond at (ap) 2 = is the value obtained after 
performing a linear extrapolation of the subtracted data. 



The renormalization factors for the one-derivative vec- 
tor and axial- vector operators, Z D v y and Z D " A , fall into 
different irreducible representations of the hypercubic 
group, depending on the choice of the external indices, 
fi, v. Hence, we distinguish between Z^y (Z^ A ) and 

Zdv (Z'da')- F° r ^ ne conversion factors from RI to MS 
we used the results of Ref. [43| for the local vector and 
axial- vector operators while for the one-derivative oper- 
ators we used the expressions of Ref. [38| . Another char- 
acteristic of these renormalization constants is that they 
depend on the renormalization scale. Thus, they need 
to be converted to the continuum MS-scheme, and for 
this we use a conversion factor computed in perturbation 
theory to 0(g A ). They are also evolved perturbatively 
to a reference scale, which is chosen to be (2 GeV) 2 . 
The results are shown in Fig. [5] both before subtracting 
the perturbative 0(a 2 )-terms and after. Using the sub- 
tracted data we find the values given in Table III. 

These are the values that we use in this work to renor- 
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Zv 
Za 

^DV 
^DV 
^DA 
DA 


0.625(2) 0.664(1) 
0.757(3) 0.771(2) 
1.019(4) 1.048(5) 
1.053(11) 1.105(4) 
1.086(3) 1.112(5) 
1.105(2) 1.119(6) 



TABLE III: Renormalization constants in the chiral limit at 
P = 1.95 and P = 2.10 in the MS-scheme at // = 2 GeV. 



malize the lattice matrix elements. The numbers in the 
parenthesis correspond to the statistical error. Our full 
results for the renormalization functions of the fermion 
field, local and one derivative bilinears along with the 
systematic error analysis will appear in a separate pub- 
lication. 



fermions. These are computed at different lattice spac- 
ings ranging from a ~ 0.1 fm to o ~ 0.06 fm. As can 
be seen, no sizable cut-off effects are observed. Lattice 
data computed using different volumes are also consistent 
down to pion masses of about 300 MeV, where we have 
different volumes. In a nutshell, our results do not indi- 
cate volume or cut-off effects larger than our current sta- 
tistical errors. A dedicated high statistics analysis using 
the Nf =2+1+1 ensemble at = 354 MeV has shown 
that contributions from excited states are negligible for 
9a [33L |34| . In recent studies, the so called summation 
method, that sums over the time-slice t where the current 
is inserted, is used as an approach that better suppresses 
excited state contributions 44]. Using this method to 
analyze lattice results at near physical pion mass it was 
demonstrated that, in fact, the value of gA decreases 0- 
These studies thus indicate that the source of the dis- 
crepancy does not seem to come from excited states con- 
tamination. 



III. LATTICE RESULTS 

In this section we present our results on the nucleon 
electromagnetic form factors, Ge(Q 2 ) and Gm(Q 2 ), and 
the axial-vector form factors, Ga(Q 2 ) and G P (Q 2 ). We 
also show the n — 2 generalized form factors for the 
one-derivative vector operator, A 2 o{Q 2 ), -B2o(Q 2 ) and 
G2o(Q 2 ), and the one-derivative axial- vector oprator, 
A20 (Q 2 ) and B2o(Q 2 )- The numerical values are given 
in the Tables in Appendix A. The dependence of these 
quantities on the momentum transfer square, Q 2 , the lat- 
tice spacing, as well as on the pion mass is examined. 
We also compare with recent results from other collabo- 
rations. 

As we already mentioned, most of the results are ob- 
tained for isovector quantities. For the renormalized nu- 
cleon matrix element of the operators with up to one 
derivative we thus consider 

<-> _ 
W7{ M D v } u - d7 {/J D V } d, 

u-y 7^ D v } u - d-y 7{ (tl D v } d , 

in the MS scheme at a scale ijl = 2 GeV. Note that the lo- 
cal vector and axial-vector operators are renormalization 
scale independent, thus the conversion to the MS scheme 
is irrelevant. 

In order to study the spin content of the nucleon we 
also compute the isoscalar matrix elements of the one- 
derivative vector operator, as well as, the isoscalar axial 
charge assuming, in all cases, that the disconnected con- 
tributions are negligible. 



A. Nucleon form factors 

In Fig. [6] we present our results for the axial charge 
gA = Ga{0) using Nf=2 and Nf =2+1+1 twisted mass 
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FIG. 6: Results for the nucleon axial charge with (i) Nf=2 
twisted mass fermions with a = 0.089 fm (filled red circles 
for L = 2.1 fm and filled blue squares for L = 2.8 fm), a = 
0.070 fm (filled green triangles), and a = 0.056 fm (open star 
for L = 2.7 fm and open square for L = 1.8 fm) [3[ (ii) 
Nf =2+1+1 twisted mass fermions with a = 0.086 fm (open 
circle) and a = 0.066 fm (square with a cross). The asterisk 
is the physical value as given in the PDG [211 ]. 



In Fig. [7] we compare our results to other recent lat- 
tice QCD data obtained with different actions. We show 
results obtained using domain wall fermions (DWF) 
clover fermions j48j, a mixed action with 2+1 flavors of 
asqtad-improved staggered sea and domain wall valence 
fermions [45[ referred to as hybrid, and 7V/=2+l of tree- 
level clover-improved Wilson fermions coupled to dou- 
ble HEX-smeared gauge fields @, S3] ■ We observe that 
all these lattice results are compatible. This agreement 
corroborates the fact that cut-off effects are negligible 
since these lattice data are obtained with different dis- 
cretized actions without being extrapolated to the con- 
tinuum limit. The recent result of Ref. 47] at almost 
physical pion mass shows about 10% deviation from the 
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FIG. 7: The nucleon axial charge for twisted mass fermions, 
Nf=2 (filled red circles) and iV/ =2+1+1 (filled blue squares), 
as well as, results using other lattice actions: Filled (green) 
triangles correspond to a mixed action with 2+1 flavors of 
staggered sea and domain wall valence fermions [iH ]. crosses 
to TV/ =2+1 domain wall fermions [5|, open triangles to Nf=2 
clover fermions [461 ] and open (cyan) circles to Nf— 2+1 of 
tree-level clover-improved Wilson fermions coupled to double 
HEX-smeared gauge fields [47t ] . 
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FIG. 8: The nucleon axial charge for twisted mass fermions 
(Nf=2 and iV/=2+l+l), as well as results using other lat- 
tice actions versus Lm^. Black symbols denote results at 
almost physical pion mass obtained using iV/=2 and 
iV/=2+l [43] clover fermions. The rest of the notation is 
the same as that in Fig. 



physical value of g e j£ p — 1.267 [2l|. This is a well-known 
puzzle and various directions have been explored to iden- 
tify the source of the discrepancy [HIM El El- In Fi gEl 
we also include the recent results obtained using Nf—2 
clover fermions at three lattice spacings a = 0.076 fm, 
0.071 fm and 0.060 fm (46[. They include a result at al- 
most physical pion mass, which is clearly higher than the 
corresponding one obtained in Ref. |47| . As already re- 
marked, the latter was shown to even decrease if one uses 
the summation method Q- In Ref. [46[ it is argued that 



volume corrections are sizable and increase the value of 
gA- We note that all lattice data shown in Fig. [7] are not 
volume corrected. In order to assess, which of these re- 
sults would suffer from large volume corrections we show 
in Fig. [5] gA as a function of Lm^. The data points 
at almost physical pion mass are shown with the black 
symbols. The result from Ref. [I?} at Lm, = 4.2 is lower 
than the one from Ref. [46[ at Lm^ — 2.74. Thus volume 
effects alone may not account for the whole discrepancy 
and therefore, there is still an open issue in the evaluation 
of gA- 
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FIG. 9: Comparison of the 7V/=2+l+l twisted mass data on 
Ga{Q 2 ) (upper) and G P (Q 2 ) (lower) for the two different pion 
masses considered. Filled blue squares correspond to /? = 2.10 
and m-r: = 210 MeV, while filled red circles correspond to 
/3 = 1.95 and = 354 MeV. The dashed lines are the dipole 
fits on the lattice data, while the solid green line is the dipole 
fit of experimental data for Ga(Q 2 ) 50] in combination with 
pion-pole dominance for G P (Q 2 ). 

Next, we study the dependence of the axial form fac- 
tors on the momentum transfer, Q 2 . In Fig.[5]wc compare 
our iV/=2+l+l results for G A (Q 2 ) and G P (Q 2 ) as the 
pion mass decreases from 354 MeV to 210 MeV. As can 
be seen, the dependence on the pion mass is very weak 
for Ga{Q 2 ) whereas for G P (Q 2 ) a stronger dependence 
is observed in particular at low Q 2 . This is not surpris- 
ing since G P (Q 2 ) is expected to have a pion-pole depen- 
dence that dominates its Q 2 -dependence as Q 2 — > 0. The 
solid line is the result of a dipole fit to the experimental 
electroproduction data for Ga(Q 2 )- Assuming pion-pole 
dominance we can deduce from the fit to the experimen- 
tal data on Ga{Q 2 ) the expected behavior for G P (Q 2 ), 
shown in Fig. As can be seen, both quantities have a 
smaller slope with respect to Q 2 than what is extracted 
from experiment. Such a behavior is common to all the 
nucleon form factors and it remains to be further inves- 
tigated if reducing even more the pion mass will resolve 
this discrepancy. The Q 2 -dependence of the lattice QCD 
data for Ga(Q 2 ) can be well parameterized by dipole 
Ansatz of the form 



G A (Q 2 



9A 



(l + Q 2 /m 2 A Y 



(22) 



10 



as it was done for the experimental results. Likewise, 
assuming pion-pole dominance we fit G P {Q 2 ) to the form 



G P (Q 2 ) 



G A (Q 2 )G p (0) 
(Q 2 +m 2 p ) 



(23) 



In both fits we take into account lattice data with Q 2 up 
to a maximum value of (1.5) 2 GeV 2 . The values of the 
parameter tua extracted from the fit for the two ensem- 
bles are 

f3 = 1.95 : m A = 1.60(5) GeV 
/3 = 2.10 : m A = 1.48(12) GeV. 

These are higher than the experimental value of m e A xp = 
1.069 GeV [50] extracted from the best dipole parameter- 
ization to the clectroproduction data. This deviation be- 
tween lattice and experimental data reflects the smaller 
slope in the lattice QCD data. Another observation is 
that the fits for G p (Q 2 ) are strongly dependent on the 
lowest values of Q 2 taken in the fit due to the strong 
Q 2 -dependence of G P (Q 2 ) at low Q 2 . 

In Figs. [10] and [TTJ we compare results using the two 
iV/=2+l+l ensembles with those obtained with Nf=2 
ensembles at similar pion masses. We do not observe 
large deviations between Nf—2 and Nf =2+1+1 result 
showing that strange and charm quark effects are small, 
as expected. 




Q- (GeV 2 ) 



FIG. 10: The Q 2 -dependence of the form factors Ga and 
G p for i) N f =2 at rrw = 377 MeV, a = 0.089 fm (filled red 
circles); ii) AT/=2+l+l at m n = 354 MeV, a = 0.086 fm. The 
solid line in the upper plot shows the resulting dipole fit to 
the experimental data on Ga(Q 2 ) [5[j. Assuming a pion-pole 
dependence for G P (Q 2 ) and using the fit on Ga{Q 2 ) shown 
in the upper panel produces the solid line shown in the lower 
panel for G P (Q 2 ). 

It is interesting to compare our TMF results to those 
obtained using different fermion discretization schemes. 
We collect recent lattice QCD results in Figs. [T2l and [T3l 
at similar pion masses. As can be seen, in the case of 
Ga{Q 2 ) there is agreement of our results with those ob- 
tained using DWF and the hybrid approach. For G P (Q 2 ) 
hybrid results obtained on a larger volume are higher at 
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FIG. 11: The Q 2 -dependence of the form factors Ga (up- 
per) and G p (lower) for Nf =2+1+1 twisted mass fermions at 
= 210 MeV, a = 0.066 fm (filled blue squares) and N f =2 
twisted mass fermions at = 262 MeV and a = 0.056 fm 
(filled red circles). The rest of the notation is the same as 
that in Fig. [T0l 



small Q 2 -values. This is an indication that volume effects 
are larger for quantities like G P (Q 2 ) for which pion cloud 
effects are expected to be particularly large at small Q 2 . 
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FIG. 12: Q 2 -dependence of Ga{Q 2 ) for V/=2+l+l at m w = 
354 MeV (filled blue squares) and the N f =2 | at m T = 
298 MeV (filled red circles) twisted mass data on a lattice 
with spatial length L = 2.8 fm and similar lattice spacing. We 
also show results with N f =2+1 DWF at m n = 329 MeV, L = 
2.7 fm (crosses) [{J and with a hybrid action with 7V/=2+l 
staggered sea and DWF at = 356 MeV and L = 3.5 fm 
(open orange circles) [45| . 



We next discuss the results obtained for the isovec- 
tor electromagnetic form factors, Ge(Q 2 ) and Gm(Q 2 )- 
In Fig. [J3 we compare our V/=2+l+l results as the 
pion mass decreases from 354 MeV to 210 MeV. As 
can be seen, the values for both quantities decrease to- 
wards the experimental values shown by the solid line, 
which is J. Kelly's parameterization to the experimen- 
tal data [5l|. In particular, for Gm(Q 2 ) lattice results 
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FIG. 13: The (^-dependence of G P (Q 2 ). The notation is the 
same as that in Fig. 1121 



at m T = 210 MeV become consistent with the experi- 
mental results. In order to extract the value of Gm(0), 
we need to extrapolate lattice results at finite Q 2 . We 
parameterized both form factors by a dipole form 



G E (Q 2 )= 



1 



G M {Q 2 



(l+Q 2 /m 2 E ) 2 ' 

Gm(0) 
\l+Q 2 lm 2 M f 



(24) 



The values of Gm (0) extracted are shown in Fig. Q3] as 
well as, the resulting fits with the dashed lines. The over- 
all trend of the lattice QCD data clearly shows that as 
the pion mass decreases they approach the experimental 
values. However, even at = 210 MeV the value of 
Gm(0), which determines the isovector anomalous mag- 
netic moment, is still underestimated. In Table IIVI we 
tabulate the resulting fit parameters m^, Gm(0) and 
tom for the two N $=2+1+1 ensembles extracted from 
the dipole fits of Eqs. (pM)) . 



p 


m E (GeV) 


G M (0) 


niM (GeV) 


1.95 


1.17(32) 


3.93(12) 


1.30(08) 


2.10 


0.86(07) 


3.86(34) 


0.99(15) 



TABLE IV: Results on the nucleon electric and magnetic mass 
extracted by fitting to the dipole form of Eq. (|24[) . 



In Fig. [15] we show the Q 2 dependence of Ge(Q 2 ) and 
G M {Q 2 ) at (3 = 2.10 and to^ = 210 MeV comparing it 
to the smallest available pion mass of 262 MeV obtained 
using Nf=2 ensembles. Once again we do not observe 
any sizable effects due to the strange and charm quarks 
in the sea. 

It is useful to compare TMF results to those obtained 
within different fcrmion discretization schemes. In par- 
ticular, we compare in Figs. [16] and I17t vith results ob- 
tained using Nf =2+1 DWF [4], Nf=2 Wilson improved 



clover fermions [48| and using the hybrid action [45[ for 
a pion mass of about 300 MeV. We see a nice agree- 
ment among all lattice results for Ge{Q 2 ), confirming 
that cut-off effects are small for these actions. In the 
case of Gm(Q 2 ) there is also an overall agreement except 
in the case of the Nf=2 clover results. These results are 
somewhat lower and are more in agreement with our re- 
sults at TOtt = 210 MeV. The reason for this is unclear 
and might be due to limited statistics as these data carry 
the largest errors. 
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FIG. 14: Comparison of the V/=2+l+l twisted mass data 
on Ge(Q 2 ) (upper) and Gm{Q 2 ) (lower) for the two different 
pion masses considered. The solid lines are Kelly's parame- 
terization of the experimental data [5l|. whereas the dashed 
lines are dipole fits to the lattice QCD data. 
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FIG. 15: The Q 2 -dependence of G E (Q 2 ) (upper) andG M (<2 2 ) 
(lower) for V/=2+l+l TMF at = 210 MeV (filled blue 
squares) and N f =2 TMF at m n = 262 MeV (filled red circles). 

Having fitted the electromagnetic form factors we can 
extract the isovector anomalous magnetic moment and 
root mean square (r.m.s.) radii. The anomalous mag- 
netic moment is given by the Pauli form factor F-zifS) 
and the slope of F\ at Q 2 = determines the transverse 
size of the hadron, (rj_) = — idFi/dQ 2 \Q2 =Q . In the non- 
rclativistic limit the r.m.s. radius is related to the slope 
of the form factor at zero momentum transfer. Therefore 
the r.m.s. radii can be obtained from the values of the 
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FIG. 16: The Q 2 -dependence of Ge{Q 2 ). We show results 
for JV/=2+l+l at = 354 MeV (filled blue squares) and 
AT/=2 | at m T = 298 MeV (filled red circles) TMF data on 
a lattice with spatial length L — 2.8 fm and similar lattice 
spacing. We also show results with Nf =2+1 DWF at = 
297 MeV, L = 2.7 fm (crosses) 0], with a hybrid action with 
AT/=2+l staggered sea and DWF at = 293 MeV and 
L — 2.5 fm (open orange circles) [45l ]. and Nf=2 clover at 
= 290 MeV and L = 3.4 fm (asterisks) The solid 

line is Kelly's parameterization of the experimental data [5ll ] 
from a number of experiments as given in Ref. [5l| . 
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FIG. 17: The Q -dependence of Gm(Q )■ The notation is the 
same as that in Fig. 1161 
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(25) 

The electric and magnetic radii are given by (r E M ) — 
12 /m E M and can be directly evaluated from the values 
of the parameters listed in Table ITVl In Fig.[TH]we present 
our results on the anomalous magnetic moment, Dirac 
and Pauli r.m.s. radii. As can be seen, the new results 
at m w = 210 MeV, although they are still lower than the 
experimental value, show an increase towards that value. 
In Ref. |2| an analysis of the results using the summation 



FIG. 18: Twisted mass fermion results with Nf=2 0] and 
with Nf =2+1+1, for the isovector anomalous magnetic mo- 
ment, K p ~ n in Bohr magnetons (upper), Dirac r.m.s. radius 
(middle) and Pauli r.m.s. radius (lower) panel. The notation 
is the same as that in Fig. [6] 



method at = 147 MeV with Nf =2+1 clover fermions 
was carried out. It was shown that the value of these 
three quantities increases to bring agreement with the 
experimental value. This is an encouraging result that 
needs to be confirmed. 
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B. Nucleon generalized form factors with one 
derivative operators 

In this section we present results on the nucleon matrix 
elements of the isovector one-derivative operators defined 
in Eq. (|22]). Like g A , A 20 (Q 2 =0) and A 20 (Q 2 =0) can be 
extracted directly from the corresponding matrix element 
at Q 2 = 0. On the other hand, B 20 (Q 2 =0), C 20 (Q 2 =0) 
and B 20 (Q 2 =0), like Gm and G p , can not be extracted 
at Q 2 = 0. Therefore one needs to extrapolate lattice 
data at Q 2 y^0 by performing a fit. 

In Fig. Ql|] we compare our lattice data of the un- 
polarized and polarized isovector moments obtained for 
Nf=2 [1] TMF for different lattice spacings and volumes 
to the N f =2+l+l TMF results of this work. As can be 
seen, there arc no detectable cut-off effects for the lattice 
spacings considered here, nor volume dependence at least 
for pion masses up to about 300 MeV where different vol- 
umes were analyzed. Also, there is consistency among re- 
sults obtained using Nf—2 and Nf =2+1+1 gauge config- 
urations indicating that strange and charm quark effects 
are small. We would like to point out that the renormal- 
ization constant for the vector one-derivative operator is 
larger by about 2% than the one used in Ref. [34| since 
in converting to MS we used the 2-loop conversion fac- 
tor instead of the 3- loop result, thus increasing the value 
of (a^u-d. As in the case of the nucleon axial charge, a 
number of studies were undertaken to examine the role of 
excited states in the extraction of (x) u -d- A high statis- 
tics analysis carried out with twisted mass fermions at 
TOtt = 354 MeV has shown that excited state contam- 
ination accounted for a decrease of about 10% in the 
value of (x) u -d as compared to the value extracted us- 
ing sink-source separation of about 1 fm [H, [Hj|. The 
most noticeable behavior regarding these TMF results 
is that the values obtained at = 210 MeV for both 
(x) u -d and (x)Au-Ad approach the physical value. We 
would like to remark that the phenomen olog ical value of 
{x)u-d extracted from different analysis [52h57| shows a 
spread, which, however, is significantly smaller than the 
discrepancy as compared to the deviation shown by lat- 
tice data for pion masses higher than the physical point. 
The same applies for (x)Au-Ad [H,[59j]. 

Recent results on A 2 q and A 2 q from a number of 
groups using different discretization schemes are shown 
in Fig. [20l We limit ourselves to results extracted from 
fitting to the ratio given in Eq. dT7|) taking a source- 
sink separation of 1 fm to 1.2 fm. Once more, there 
is an overall agreement among these lattice data indi- 
cating that cut-off effects are small for lattice spacings 
~ O.lfm, for the improved actions used. The decrease 
seen using TMF at = 210 MeV is corroborated 
by other recent results at near physical pion masses: 
for (x) u -d results from Ref. [6] using clover- improved 
fermions at = 157 MeV and Lm, = 2.74, as well 
as, from Ref. using Nf =2+1 flavors of tree- level 
clover- improved Wilson fermions coupled to double HEX- 
smeared gauge fields at = 149 MeV and Lm^ = 4.2, 
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FIG. 19: Results for {x) u -d (upper) and (x)Au-Ad (lower) 
using Nf=2 and Nf =2+1+1 twisted mass fermions as a 
function of the pion mass. We show results for (i) Nf=2 
twisted mass fermions with a = 0.089 fm (filled red circles 
for L = 2.1 fm and filled blue squares for L = 2.8 fm), 
a = 0.070 fm (filled green triangles), and a = 0.056 fm (open 
stars for L = 2.7 fm and open square for L = 1.8 fm); (ii) 
Nf =2+1+1 twisted mass fermions with a = 0.0863 fm (open 
circle) and a = 0.0657 fm (square with a cross) . The physical 
point, shown by the asterisk, is from Ref. [55|] for the unpo- 
larized and from Ref. [5S| for the polarized first moment. 



also decrease towards the physical value. Furthermore, 
for the latter case, three sink-source separations up to 
1.4 fm were utilized to apply the summation method re- 
ducing the value shown in Fig. [20] further to bring it into 
agreement with the experimental one Note that this 
is opposite to what was found for g A where its value de- 
creased further away from the experimental value. The 
agreement between the values found in Refs. @ and [47j | . 
despite the different volumes, indicates that the volume 
dependence of this quantit y is small, again different from 
what was claimed in Ref. [46[ for g A - In Ref. [47| it was 
demonstrated that contributions from excited states in- 
crease as the pion mass decreases towards its physical 
value indicating that excited state contamination may 
explain the discrepancy between lattice results and the 
experimental value. Further studies of excited state con- 
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FIG. 20: Results for (x) u -d (upper) and (x)Au-Ad (lower) 
obtained in this work are shown with the red filled circles 
for Nf—2 and with the blue filled squares for Nf =2+1+1. 
We compare with (i) iV/=2+l DWF for a = 0.114 fm [g3|; 
(ii) Nf =2+1 using DWF for the valence quarks on stag- 
gered sea [45[ with a = 0.124 fm; (iii) Nf=2 clover with 
a — 0.075 fm |6l| . For (x) u -d we also show recent results 
using N f =2 clover with a = 0.071 fm 6] and Nf=2+1 of 
tree-level clover-improved Wilson fermions coupled to double 
HEX-smeared gauge fields with a = 0.116 fm 



tamination at near physical pion mass will be essential 
in order to establish this conclusion. 

The Q 2 -dependence f A2o(Q 2 ) and ^2o(Q 2 ) is shown 
in Fig.[5T]for our two Nf— 2+1+1 ensembles and for the 
Nf—2 ensemble with the smallest available mass, namely 
262 MeV. Since strange and charm quark effects have 
been shown to be small, one can study the dependence on 
the pion mass by comparing with results obtained using 
Nf—2 TMF. As the pion mass decreases from 354 MeV 
to 262 MeV there is no significant change in the values of 
A2o{Q 2 ) and A 2 o(Q 2 ) over the whole Q 2 range. Reduc- 
ing the pion mass further to 210 MeV leads to a larger 
decrease in the values of both A2o(Q 2 ) and A20 (Q 2 ) in- 
dicating that near the physical regime the pion mass de- 
pendence becomes stronger. Such a pion mass depen- 
dence is what one would expect if the lattice QCD data 
at Q 2 = are to agree with the experimental value. In 
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FIG. 21: The Q 2 -dependence of A 20 {Q 2 ) (upper) and A 20 {Q 2 ) 
(lower) for TV/ =2 with a = 0.056 fm and m n = 262 MeV (filled 
green diamonds), and Nf =2+1+1 with i) a = 0.066 fm and 
= 210 MeV (filled blue squares); ii) a = 0.086 fm and 
tutt — 354 MeV (filled red circles). 
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FIG. 22: The Q 2 -dependence of A 20 (Q 2 ) (upper) and A 20 (Q 2 ) 
(lower) shown for i) Nf—2 twisted mass fermions for a = 
0.089 fm, m w = 377 MeV (filled red circles) 4]; ii) JV/=2+l+l 
twisted mass fermions (this work) for a — 0.086 fm and = 
354 MeV; iii) Nf—2 clover fermions for a ~ 0.08 fm and ~ 
350 MeV (open cyan diamonds) [62]|: and iv) Nf =2+1 with 
DWF valence on a staggered sea for a — 0.124 fm and m n = 
356 MeV (open orange circles) [451 ] . 



Fig. [22] we compare our results using TMF to hybrid re- 
sults and, for A2o(Q 2 ), we also include Nf=2 clover at 
similar pion masses. There is an overall agreement be- 
tween clover and TMF for A2o(Q 2 ), whereas the hybrid 
data are somewhat lower. The fact that they are renor- 
malized perturbatively might explain their lower values. 
Before closing this section we present in Fig. |2"31 results 
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FIG. 23: The Q 2 -dependence of B 20 (Q 2 ), C 20 (Q 2 ) and 
B 20 (Q 2 ) for N f =2+1+1 computed at /? = 1.95 (ra* = 354 
MeV) and /3 = 2.10 (m, = 210) MeV. The dashed lines show 
the linear fits to B2o(Q 2 ), C 2 o(Q 2 ) and B 2 o(Q 2 ) to extract 
the value at Q 2 — shown here. 



for S 20 (Q 2 ), C 20 {Q 2 ), B 20 (Q 2 ) for the two N f =2+1+1 
ensembles. All these three GFFs can not be extracted 
at Q 2 =0 directly from the matrix element and therefore 
we must extrapolate them using an Ansatz to fit the in- 
dependence. We performed two types of fits: a linear and 
a dipole fit. Note that for small Q 2 the two are equiv- 
alent. It was generally found that a linear fit describes 
well the data with smaller errors on the fit parameters. 
We therefore use the fitted values extracted from the lin- 
ear fit summarized in Table fVl C2o(Q 2 ) is consistent with 
zero for all values of Q 2 . 



p 


B 20 (0) (GeV) 


C 20 (0) 


B 2 o(0) 


1.95 


0.344(19) 


-0.009(09) 


0.648(71) 


2.10 


0.205(62) 


0.016(34) 


0.518(251) 



TABLE V: Results on B 20 (Q 2 = 0), C 20 (Q 2 
B 2 q(Q 2 = 0) by fitting to a linear Q 2 -dependence. 
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IV. PROTON SPIN 

How much of the proton spin is carried by the quarks 
is a question that is under study ever since the results of 
the European Muon Collaboration (EMC) claimed that 
the quarks carried only a small fraction of the proton 
spin [63j . This became known as the "proton spin crisis" . 
It was proposed that gluons in a polarized proton would 
carry a fraction of the spin, which however would be un- 
naturally large if it were to resolve the EMC spin crisis. 
It is now understood that the resolution of this puzzle 



requires to take into account the non-perturbative struc- 
ture of the proton • In order to use our lattice results 
to obtain information on the spin content of the nucleon 
we need to evaluate, besides the isovector moments, the 
isoscalar moments A^ d and B^ d since the total angular 
momentum of a quark in the nucleon is given by 



j q = o (^!o(o) + syo)) • 



(26) 



As already discussed, the total angular momentum J 9 
can be further decomposed into its orbital angular mo- 
mentum L q and its spin component AS 9 as 

1 



-A£ 9 + L q . 



(27) 



The spin carried by the u- and d- quarks is determined 
using AE" +d = A"q , and therefore we need the isoscalar 
axial charge. The isoscalar quantities take contribu- 
tions from the disconnected diagram, which are notori- 
ously difficult to calculate and are neglected in most cur- 
rent evaluations of GFFs. These contributions are cur- 
rently bein g co mputed using improved stochastic tech- 
niques [65l . |66| . Under the assumption that these are 
small we may extract information on the fraction of the 
nucleon spin carried by quarks. 

In Fig. [24] we show our results for the isoscalar 
G A (Q 2 r+ d , A 20 (Q 2 r+ d , B 20 (Q 2 r +d and C 20 (Q 2 ) u+d 
for the two Nf=2+1+1 ensembles analyzed in this work. 
It was shown using the Nf=2 ensembles at three lattice 
spacings smaller than 0.1 fm [l| that cut-off effects are 
small. We expect a similar behaviour for our A/=2+l+l 
ensembles. Therefore, we perform a chiral extrapola- 
tion using directly all our lattice data for the Nf=2 and 
Nf =2+1+1 ensembles. Having both isoscalar and isovec- 
tor quantities we can extract the angular momentum J u 
and J d carried by the u- and d- quarks. In order to ex- 
tract these quantities we need to know the value of B 20 
at Q 2 = 0. As explained already, one has to extrapolate 
the lattice results using an Ansatz for the Q 2 -dependence 
to extract B 2 o at Q 2 = and two ansatze were consid- 
ered for the Q 2 -dependence, a dipole and a linear form. 
For the linear fit we use two fitting ranges one up to 
Q 2 = 0.25 GeV 2 and the other up to Q 2 = 4 GeV 2 . 
Thus the extrapolation introduces model dependence in 
the extraction of the quark spin J q . The values of B 2 q 
extracted using these three ansatze are consistent, with 
the dipole fit resulting in parameters that carry large er- 
rors. In extracting the angular momentum we thus use 
the data extracted using the extended range linear fit and 
given in Table [V] 

We first compare in Fig. [25] our results for the u- and 
d- quark angular momentum J 9 , spin AS 9 and orbital 
angular momentum L q to those obtained using the hy- 
brid action of Ref. j45] . As can be seen, the lattice data 
are in agreement within our statistical errors indicating 
that lattice artifacts are smaller than the current statis- 
tical errors, also for these quantities. In order to get 
an approximate value for these observables at the phys- 
ical point we perform a chiral extrapolation using heavy 
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FIG. 24: The Q 2 -dependence of the isoscalar Ga(Q 2 ), 
A 20 {Q 2 ) and B 20 {Q 2 ) for N f =2+1+1 computed at /3 = 1.95 
{m n = 354 MeV) and /3 = 2.10 (m n = 210) MeV. 



baryon chiral perturbation theory (HB^PT). Combining 
the expressions for A20 and B20 [6jg, |69( in the isoscalar 
and isovector cases we obtain the following form for the 
angular momentum 



jq _ q "-7T 1 , q 2 



(28) 



and take A 2 = 1 GeV 2 . We also carry out a chiral fit 
using Q(p ) covariant baryon chiral perturbation theory 
(CBxPT) 70]. All the expressions are collected in Ap- 
pendix B for completeness. As noted these chiral extrap- 
olations are to give an indicative idea of what one might 
obtain since their range of validity may require using pion 
masses closer to the physical point. 

In order to correctly estimate the errors both on the 
data points and on the error bands, we apply an ex- 
tended version of the standard jackknife error procedure 
known as super jackknife analysis [45| . This generalized 
method is applicable for analyzing data computed on sev- 
eral gauge ensembles. This is needed for carrying out the 
chiral extrapolations for the angular momentum J q , or- 
bital angular momentum L q and spin AS 9 . Although, 
there is no correlation among data sets from different 
gauge ensembles, the data, within each ensemble arc cor- 
related. This analysis method allows us to consider a 
different number of lattice QCD measurements for each 
ensemble taking into account correlations within each en- 
semble correctly. It should be apparent that the super- 
jackknife reduces to the standard jackknife analysis in 
the case of a single ensemble. 

In Fig. [26] we show the chiral fits for J q . In the upper 
panel we show the chiral extrapolation using CB%PT and 
in the lower the extrapolation using HBxPT. Both have 
the same qualitative behavior yielding a much smaller 
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FIG. 25: Comparison of TMF results (filled symbols) to those 
using a hybrid action [451 ] (open symbols). The upper panel 
shows the angular momentum J u and J for u- and d- quarks 
respectively (blue filled squares for A^/=2+l+l and filled red 
circles for Nf—2). The lower panel shows the quark spin 
(same symbols as for J q ) and the orbital angular momentum 
(filled green triangles for Nf=2 and filled magenta diamonds 
for Nf— 2+1+1). The errors are determined by carrying out a 
superjacknife analysis described in Ref. [4^]. The experimen- 
tal value of AT, u ' d is shown by the asterisks and are taken 
from the HERMES 2007 analysis [13] ■ 



contribution to the angular momentum from the d-quark 
than that from the u-quarks. In the plot we also show 
the band of allowed values if the fit were performed on 
data that used the Q 2 = extrapolated values of B20 
from the limited range linear fit. As can be seen, the 
two bands are consistent. Had we used a dipole Ansatz 
for the Q 2 = extrapolation, the error band would also 
be consistent but much larger, especially for smaller pion 
masses, where there are no lattice data. Therefore, for 
the rest of the discussion we only show the extrapolation 
bands obtained using the limited and full Q 2 range linear 
fits. These results are in qualitative agreement with the 
chiral extrapolations using the data obtained with the 
hybrid action [iBj . 
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FIG. 26: Chiral extrapolation using CB^PT (upper) and 
HB^PT (lower) for the angular momentum carried by the 
u-and d- quarks. The red band is the chiral fit using the data 
for i?2o(Q 2 = 0) obtained by a linear extrapolation of B2o(Q 2 ) 
using Q 2 values up to Q 2 = 4 GeV 2 whereas the green band 
is the fit using values of £?2o(0) extracted from a linear ex- 
trapolation of i?2o(Q 2 ) using Q 2 values up to ~ 0.25 GeV 2 . 
The data shown in the plot are obtained from the extended 
linear Q 2 extrapolation. Filled red circles are data for Nf—2 
at /3 = 3.9, filled green triangles for N f =2 at f) = 4.05, filled 
magenta diamonds for Nf—2 at /3 = 4.2, filled light blue in- 
verted triangle for TV/ =2+ 1+1 at p = 1.95 and filled blue 
square for iV/=2+l+l at p = 2.10. 



In Fig. [57] we show separately the orbital angular mo- 
mentum and spin carried by the u- and d- quarks. The 
total orbital angular momentum carried by the quarks 
tends to small negative values as we approach the phys- 
ical point. This is a crucial result and it would be 
important to perform a calculation at lower pion mass 
to confirm that this trend towards negative values re- 
mains [71| . After chiral extrapolation, the value obtained 
at the physical point is consistent with zero in agreement 
with the result by LHPC. We summarize the values for 
the angular momentum, orbital angular momentum and 
spin in the proton at the physical point in Table fVTl The 
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FIG. 27: Chiral extrapolation using HBxPT. The upper 
graph shows the spin and orbital angular momentum carried 
by u- and d- quarks, whereas the middle and lower graphs 
show the spin and orbital angular momentum carried sepa- 
rately by the u- and d- quarks. The errors are determined 
through a superjacknife analysis. The physical points, shown 
by the asterisks are from the HERMES 2007 analysis [g^]. 
The notation is the same as that in Fig. 1261 
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T[T) , |)| 

rirSx^ 1 


experiment 


— J u d — 




0.218179(12) 




TU-\-d 
J 




0.168(42) 




J U 


0.249616(7) 


0.173(35) 




rd 
J 


-0.000242(6) 


-0.003(17) 




A \^u — d /r> 
ZAZj 1 1 




0.462(9 j 












AE"/2 




0.440(29) 


0.421(6) 


AE d /2 




-0.137(19) 


-0.214(6) 


L u-d 




-0.252(8) 








-0.141(48) 




L u 




-0.196(23) 




L d 




0.055(24) 





TABLE VI: Values of nucleon spin observables at the ph ysic al 
point using CBxPT and HBxPT and from experiment [671 ] - 



value obtained for the AS" agrees with the experimental 
value whereas the extrapolated LQCD value of AE d is 
less negative. As already pointed out, results closer to 
the physical pion mass will be essential to resolve such 
discrepancies. In addition, the computation of the dis- 
connected diagrams will eliminate a remaining system- 
atic error and will enable us to have final results on the 
spin carried by the quarks and consequently on the gluon 
contribution to the nucleon spin. 

V. CONCLUSIONS 

We have performed an analysis on the generalized form 
factors G E (Q 2 ), G M (Q 2 ), G A (Q 2 ), G P (Q 2 ), A 20 (Q 2 ), 
B 20 (Q 2 ), C 2a {Q 2 ), A 2a (Q 2 ) and B 20 (Q 2 ), extracted 
from the nucleon matrix elements of the local and 
one-derivative vector and axial-vector operators using 
Nt=2+1+1 flavors of twisted mass fermions. Our re- 
sults are non-perturbatively renormalized and they are 
presented in the MS scheme at a scale of 2 GeV. The 
comparison of the results using Nf=2 and Nf— 2+1+1 
twisted mass fermions with the results obtained using 
other discretizations show an overall agreement for pion 
masses down to about 200 MeV. The compatibility of 
Nf~2 data with those including a dynamical strange 
and a charm quark is an indication that any systematic 
effect of strange and charm sea quark effects on these 
quantities for which disconnected contributions were ne- 



glected is small. The twisted mass fermion results on 
the axial nucleon charge remain smaller than the exper- 
imental value. The recent results using Nf=2 [46[ and 
Nf =2+1 0] clover- improved fermions near the physical 
pion mass are somewhat in conflict with each other and 
hard to interpret in a consistent way. Therefore, fur- 
ther investigation is required to resolve the issue. For 
the unpolarized isovector momentum fraction lattice re- 
sults show a decrease as we approach the physical pion 
mass with indications of excited state contamination that 
needs further investigation. 

We also analyze the corresponding isoscalar quantities 
using directly our lattice data. Of particular interest here 
is to extract results that shed light on the spin content of 
the nucleon. Assuming that the disconnected contribu- 
tions to the isoscalar quantities are small we can extract 
the spin carried by the quarks in the nucleon. For the chi- 
ral extrapolations of these quantities we use HBxPT and 
CBxPT theory applied to all our N f =2 and AT) =2+1+1 
data. We find that the spin carried by the d-quark is al- 
most zero whereas the u-quarks carry about 50% of the 
nucleon spin. This result is consistent with other lattice 
calculations [45j . 
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Appendix A: Numerical results for the isovector sector 



(GeV) 


(Q) 2 


Ge 


Gm 


Ga 




(no. confs) 












P = 1.95, 32 3 x 64 




0.0 


1.000(1) 


3.930(117) 


1.141(18) 


18.211(9.209) 




0.192 


0.734(6) 


2.979(61) 


0.995(14) 


9.462(399) 


0.354 


0.372 


0.570(7) 


2.355(46) 


0.872(12) 


6.116(226) 


(950) 


0.542(1) 


0.469(10) 


1.937(47) 


0.775(14) 


4.512(209) 




0.704(1) 


0.392(12) 


1.676(57) 


0.714(21) 


3.117(208) 




0.859(2) 


0.331(11) 


1.405(46) 


0.642(18) 


2.591(134) 




1.007(2) 


0.288(13) 


1.250(53) 


0.589(21) 


2.134(129) 




1.287(3) 


0.208(20) 


0.950(79) 


0.480(39) 


1.441(182) 




1.420(4) 


0.185(20) 


0.865(85) 


0.450(41) 


1.249(163) 




= 2.10 


48 3 x 96 










0.0 


1.006(6) 


3.855(342) 


1.164(62) 


14.880(11.790) 




0.147 


0.722(21) 


2.849(198) 


1.034(47) 


10.454(1.445) 


0.210 


0.284 


0.565(23) 


2.347(142) 


0.909(42) 


6.317(783) 


(900) 


0.414(1) 


0.430(30) 


1.950(153) 


0.850(52) 


5.227(699) 




0.537(1) 


0.444(41) 


1.622(170) 


0.690(68) 


2.466(723) 




0.655(2) 


0.318(29) 


1.338(120) 


0.689(53) 


2.628(395) 




0.768(3) 


0.266(32) 


1.291(136) 


0.707(71) 


2.763(481) 




0.980(4) 


0.218(52) 


1.104(237) 


0.558(129) 


2.466(701) 




1.081(5) 


0.186(44) 


0.686(164) 


0.437(110) 


1.714(541) 



TABLE VII: Results on Ge, Gm-, Ga and G p form factors at /3 = 1.95 (32 3 x 64) and /3 = 2.10 (48 3 x 96). 
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m-n (GeV) 


(Q) 2 


A 20 


B20 


A 20 


B20 


(no. confs) 












/3 = 1.95, 32 3 x 64 




0.0 


0.270(5) 


0.344(19) 


0.302(5) 


0.648(71) 




0.192 


0.242(4) 


0.292(15) 


0.281(5) 


0.582(121) 


0.354 


0.372 


0.222(4) 


0.266(15) 


0.264(5) 


0.578(65) 


(950) 


0.542(1) 


0.207(5) 


0.266(15) 


0.249(6) 


0.495(73) 




0.704(1) 


0.195(6) 


0.213(20) 


0.231(7) 


0.236(81) 




0.859(2) 


0.177(6) 


0.209(15) 


0.219(7) 


0.319(46) 




1.007(2) 


0.163(8) 


0.192(16) 


0.202(9) 


0.294(47) 




1.287(3) 


0.152(14) 


0.169(25) 


0.170(15) 


0.200(63) 




1.420(4) 


0.134(14) 


0.134(21) 


0.164(16) 


0.165(53) 


/3 = 2.10, 48 3 x 96 




0.0 


0.228(18) 


0.205(62) 


0.251(19) 


0.518(251) 




0.147 


0.206(14) 


0.184(63) 


0.242(14) 


0.793(475) 


0.210 


0.284 


0.190(12) 


0.233(52) 


0.247(14) 


0.830(244) 


(900) 


0.414(1) 


0.165(16) 


0.224(58) 


0.229(18) 


0.526(259) 




0.537(1) 


0.176(23) 


0.159(63) 


0.204(25) 


-0.446(289) 




0.655(2) 


0.152(16) 


0.159(49) 


0.180(19) 


-0.036(145) 




0.768(3) 


0.167(20) 


0.205(53) 


0.181(24) 


0.101(160) 




0.980(4) 


0.173(40) 


0.305(90) 


0.164(43) 


0.371(233) 




1.081(5) 


0.142(33) 


0.164(65) 


0.122(35) 


0.106(166) 



TABLE VIII: Results on A 20 , B 20 , A 20 and B20 form factors at j3 = 1.95 (32 3 x 64) and /3 = 2.10 (48 3 x 96). 



Appendix B: Expressions for the chiral extrapolation of the quark spin and angular momentum 



In this Appendix we collect the expression used to extrapolate our lattice data for the quark spin to the physical 
point. Throughout, we use A 2 = 1 GeV 2 , f v = 0.0924 GeV and g A = 1.267. 

In HBxPT the expressions for A 2 o(0) and i?2o(0) for the isovector combination are given by 
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and for the isoscalar by 



^o =o (o) = i?r (o) 
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The spin carried by the quarks is given by the axial coupling gA or Ai O (0) as 
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The corresponding expressions for A w (0) in the isoscalar and isovector cases are 
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For the total spin J we have 
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and the expression for AS 9 , L q and J q are of the form 



or 
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(40) 



where Q = J, AS, L. 

We also use covariant baryon chiral perturbation theory (CB^PT) for ^2o(0), -820(0), 620(0) in the isovector case 
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and the isoscalar case 
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We then extract the total spins using 
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